function [hq,eq,jq] = calcAhVec(invP,etaModel,hv,ev,jv)

k = 1/invP.disT(1,2);

%Mass matricies
T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - etaModel)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - etaModel)));
M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

hq = cell(invP.nt,1);
eq = cell(invP.nt,1);
jq = cell(invP.nt,1);

hq{1} = M*hv{1} + invP.mesh.CURL'*ev{1};
eq{1} = invP.mesh.CURL*hv{1} - jv{1};
jq{1} = S*ev{1} - T*jv{1};

for ii = 2:invP.nt
    hq{ii} = M*hv{ii} + invP.mesh.CURL'*ev{ii} - M*hv{ii-1};
    eq{ii} = invP.mesh.CURL*hv{ii} - jv{ii};
    jq{ii} = S*ev{ii} - T*jv{ii} - St*ev{ii-1} + Tt*jv{ii-1};
end